a_initial = [0.1057, 0.1420, 0.1664, 0.6209, 0.5737, 0.0520, 0.9312, 0.7286, 0.7378, 0.0634, 0.8604, 0.9344, 0.9843, 0.8589, 0.7855, 0.5133;
    0.1776, 0.3985, 0.1339, 0.0308, 0.9391, 0.3013, 0.2955, 0.3329, 0.4670, 0.6481, 0.0252, 0.8422, 0.5590, 0.8540, 0.3478, 0.4460];

% 假设的最优权重
a_d=ones(size(a_initial));

% 计算初始权重误差
e_a_initial = a_initial - a_d;

% 初始协方差矩阵 P（这里我们简化为对角矩阵，对角线上的值是权重误差的平方）
P_initial = diag(e_a_initial(:) .^ 2);

R = 1;% 观测噪声协方差
xi_k=ones(size(a_initial)).*0.01;

% 初始状态
a_k = a_initial;
P_k = P_initial;

% 模拟时间步长
num_steps = 100;
dt = 0.01;

% 循环模拟
for k = 1:num_steps
    % 假设的观测模型雅可比矩阵 O_k（这里简化为单位矩阵）
    O_k = eye(size(P_k));

     % 卡尔曼增益
      % 计算卡尔曼增益
    K_k=P_k*O_k'/(R+O_k*P_k*O_k');

     % 假设的观测模型
    o_k = a_k + xi_k;

     % 权重更新
    a_k = a_k + K_k * (o_k - O_k * a_k);

    % 协方差矩阵更新
    P_k = P_k - K_k * O_k * P_k;

     % 计算权重误差
    e_a = a_k - a_d;

    % 计算Lyapunov函数变化量
    delta_E_k = -e_a' * P_k^(-1) * e_a + r_k + (3/4) * (xi_k^2);
    
    % 更新Lyapunov函数
    E_k_plus_1 = E_k + delta_E_k;

     % 更新Lyapunov函数值
    E_k = E_k_plus_1;
end